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Two computer-aided algorithms for the design of all-pass digital filters 
are presented. The first technique is based on a linear programming ap- 
proach to solving the approximation problem posed by the minimax design 
of an all-pass digital filter. A new iterative algorithm with stability con- 
straints is offered for direct form design. The second technique implements 
a gradient search for those quadratic factors of an all-pass transfer function 
that lead to a locally optimal approximation (in the least-squares sense) 
of a desired phase function. New initial guess procedures and the parame- 
terization of linear-phase offset enhance the least-squares design procedure. 
Examples illustrating the result of both procedures are included. 

I. INTRODUCTION 

The increasing availability of digital signal processors such as those 
described in Refs. 1 and 2 has generated much interest in the algo- 
rithmic design of digital filters. One particular class of recursive 
digital filters commonly referred to as all-pass digital networks has 
an important and interesting design problem associated with it. That 
is, the design objective for this type of filter involves the following 
transfer function 

N 

H(r*) = k -^ • 

L b k z- k 

Because of the relationship between numerator and denominator 
polynomials, the number of degrees of freedom in filter design has 
been reduced to N from the usual 2N. Since the magnitude function 
of H(z~ l ) is precisely 1.0 on the unit circle, the design problem is 
focused directly on the phase variation of H(z~ l ). The importance of 
this design problem does not arise from an academic viewpoint. 

There are signal processing applications in which an influential 
factor in signal fidelity is the amount of phase distortion present in 
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Fig. 1 — (a) Original pulse, (b) Phase-distorted pulse. 

the medium. The effects of phase distortion in communication systems 
are illustrated in Refs. 3 and 4. Apart from nonlinear phase equalization 
applications, all-pass networks can be used to provide a constant 
phase shift over a specified frequency band or bands. The Hilbert trans- 
former commonly found in bandpass modulation systems is just one 
example of this application. In constructing phased arrays in radar 
and seismic research, constant phase shifters are also found to be 
useful. 6,6 Figure 1 illustrates how a constant phase offset can shape 
(or distort) the impulse response of a system where f(t) and f*(t) 
differ by a constant phase offset of r/2. A constant phase shift of any 
amount besides an odd multiple of r/2 will produce a pulse with a 
single large lobe. Equalization of this type of distortion is again 
possible by all-pass networks. 

Previous work 7 has addressed the envelope delay design problem. 
In many cases, this is sufficient but, as seen above, there are applica- 
tions where the phase function must be treated directly. 

Our design techniques are for all-pass structures where the design 
criteria stem from the phase function directly. The first technique, 
described in Section II, is a new method for designing all-pass networks 
using linear programming. This approach allows for fast (at least 
quadratic), always convergent design of phase networks. For the 
first time, stability can be treated directly in the design procedure. 
The second algorithm is based on a gradient search procedure on a 
least-squares criterion. The basic approach is analogous to those 
described in Refs. 7 to 9. The all-pass structure reduces the number 
of variables and simplifies the gradient calculations. In addition to 
developing the algorithm, we provide initial guess procedures and 
linear-phase offset parameters that enhance the algorithm. These 
initial guess procedures are new noniterative filter designs that can 
serve as excellent all-pass approximations in their own right. 

II. A LINEAR PROGRAMMING APPROACH 

A need for fast, reliable design of all-pass digital filters has been 
shown in the previous section. Linear programming techniques have 
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been found to be useful in rational function approximations 1011 and 
have been applied to the magnitude-only design of digital filters. 1213 
Here we show how the all-pass structure in digital filters can be trans- 
formed into a problem that can be handled by linear programming 
techniques also. As we shall see, the rational function differs from the 
magnitude-only case. Most importantly, this technique allows the 
question of stability to be handled directly in the design procedure. 
Other techniques that consider the phase or envelope delay variation 
of the digital filter (see Refs. 7 and 9 and Section III of this paper) 
deal with stability with a more heuristic approach. 

To develop the linear programming design method, we first recall 
that the all-pass transfer function is 

„ f _ n Pjz- 1 ) _ b N + Zfr-iz- 1 + &AT-2Z- 2 + • ■ • + boZ- N 
{Z } " QGr*) " " b Q + biz" 1 + • • • + b N z~» (la; 



P(z-i) z-»(b N z N + bar-iz"- 1 + • • • + 60) 



(lb) 



Q(z~ l ) (b N z~ N + b N . lZ ~ N+1 + • • • + b ) 
Hence, the phase function of (1) on the unit circle is 

(z N ££3 ) I = - 2*[Q( 2 - 1 )] I N=1 . (2) 

From (2) we note that the phase variation of H(z~ l ) is equivalent 
(modulo a constant multiplier and an TV sample delay term) to the 
phase of Q{z~ l ). Henceforth, we address the problem of synthesizing 
Q{z~ l ). The phase variation of Q{z~ l ) on the unit circle is 

N 

— £ bk sin 2irkf 

0[<2(2- 1 )]Im-i = tan- 1 — ^ 

£ b k cos2Trkf 
t=o 

tan 0[<2(2T»)]| 1.1-1 = Imag [Q(e-> 2 ^)]/Real [g(*-*9 J (3) 

Further, 

N 

- £ b k sin2irkf 

tan *[«(*-*")] = -^~ — = ftn' (4) 

£ 6* cos 2dfc/ UJ 

k=0 



Our design criterion is chosen to be 

Rifn) 



mm max 

(6*) n 



D(fn) ~ 



n = 0, 1, 2, •••, M, 



S(fn) 

where D(f) is the tangent desired phase function and M is a number 
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of frequency points* 0$>N) chosen to ensure adequate approximation 
over a subinterval of |/| ^ §, namely, < / < /i- • • < Jm < £• We 
recall that the desired phase function has been scaled down by — \ 
because of the factor of — 2 appearing in (2) and will have a delay of N 
samples inherent in its design by the z~ N factor of (lb). It is important 
to note here that the norm is applied to the tangent of the desired 
phase function instead of the desired phase function itself.* 

If we prevent S(f) from assuming the value zero, we seek the 
minimum value of A, 

\D(f n )S(f n ) - R(f n )\ ^AS(f n ). (5) 

Using the differential correction idea of Ref. 10, we expand the right- 
hand side of (5) in an iterative form : 

AS(fn) « AkS k (f n ) + (A - A k )S k (f n ) + [S(/ n ) - £*(/„)]A fc . (6) 

The intention is to iterate toward those values of [bj] that minimize A. 
The subscript k indicates the fcth iteration. We then have, from (5) 
and (6), 

\D(f n )S(fn) ~ R(fn)\ ~ A k S(f n ) - (A - A k )S k (f n ) =g 0, 

which translates into a familiar pair of equations 10 

[D(/ n ) + A k -]S(f n ) - R(f n ) + (A - A k )S k (f n ) ^ (7) 

[-£>(/„) + A k lS(f n ) + R(f n ) + (A - A k )S k (f n ) £ 0. (8) 
Substituting the series forms for R (f n ) and S(f n ), we have 

L {[#(/n) + A*]cos27rif n + sin2TJfn}b j + (A - A k )S k (f n ) 

^ - D(f n ) - A k (9) 

E {[-#(/») + A fc ] cos 2iri/-„ - sin2T7/„}& i + (A - A k )S k (f n ) 

^D(f n )-A k , (10) 

where 6 — 1 is the normalization made. We have in (9) and (10) an 
over-determined set of 2M equations in iV + 1 variables. The objective 
is to minimize A, one of the variables. It would seem that the condition 
S(fn) > would be necessary to solve (9) and (10). But the phe- 



* An extension into a weighted criterion can be handled, but is suppressed in this 
presentation. M was chosen to be in the range 4JV (N large) ^ M ^ ION (N small) 
in our implementation of the algorithm. 

* Therefore, the nonlinear nature of the tangent transformation may inhibit 
designs in the neighborhood of ir. 
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nomenon experienced in Ref. 10 occurs here also. That is, if £*(/„) > 0, 
^ n ^ M , then «S*+i(/») > 0, ^ n ^ M, also. 

Standard linear programming techniques can now be used on (9) 
and (10) to iterate toward a minimum A. However, no restriction has 
been made on the locations of the zeros of Q(z~ l ). Now there exist 
sufficient conditions for stability that can be written as linear con- 
straints. We have looked at two of these, e.g., a restriction that 
&ij &2, • • • , b N of (1) form a monotonic sequence 14 or the restriction that 
the sum ££=<, bit cos 2*kf > 0, V/E[0, O" (The formulation of the 
linear programming problem gives us this condition on the subset 
of [0, |] over which we are approximating.) For an example of a filter 
designed using this technique and the latter constraint to assure 
stability, refer to Fig. 2. Curve B is the sixth-order approximation to 
Curve A (only approximated over [0.075, 0.425]). 

However, the filter designer may decide that these types of con- 
straints are too restrictive for his particular applications. Nonlinear 
stability constraints, such as those found in Ref. 14, Chapter 3, can 
be included via the cutting planes algorithm, 16 but this may require 
excessive computation times. Another suggestion involves interrupting 
the standard simplex method for solving the linear programming 
problem after each iteration. We may then further constrain the b 
vector used in the next basic feasible solution to a choice (i.e., some 




FRACTION OF THE NYOUIST RATE 
NORMALIZED FREQUENCY 



Fig. 2 — Sixth-order approximation using linear programming method. 
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NORMALIZED FREQUENCY 

FRACTION OF THE NYQUIST RATE 



0.015 



0.200 




Fig. 3 — Phase error vs bandwidth for various orders of Hilbert transformers. 

"maximum") from among those vectors that would result in a stable 
filter in addition to the normal improvement of an object function. 

Using the standard formulation of the problem with no additional 
constraints or techniques necessary to assure stability, we were able 
to design many Hilbert transformer filters.* Figure 3 shows the relation- 
ship between the maximum error (recall that the tangent of the 
desired function is approximated) and a bandwidth (the filters were 



* fir designs of Hilbert transformers are well documented (see Ref. 17). There, 
90-degree phase is guaranteed, and the magnitude of 1.0 is approximated. 
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Fig. 4 — Tenth-order Hilbert transformer. 

designed* over [/, 0.5 -/],/= 0.075, 0.05, • • •, 0.225) for various 
orders of filters N = 4, 6, 8, 10. The log of the maximum error is given 
in the figure. 

The minimax approximation formulated here is performed on the 
tangent of the desired phase and not on the desired phase itself. For 
very good approximations, however, no penalties seem apparent. We 
have briefly looked at methods to design minimax phase approxima- 
tions based on the algorithm we have presented here. Our conclusions 
are that a two-stage design algorithm is required to iteratively locate 
a proper weight function that will "prewarp" the "tangent" design 
so that the weighted "tangent" design is minimax and the phase 
approximation is itself equiripple. 

Figure 4 illustrates the effect of the tangent transformation in this 
design procedure. In this figure, we see the phase of the resultant 
design (and its error function) . This is a 10th order approximation to a 
90-degree phase shift over [0.05, 0.45]. While the design guarantees 
a minimax solution (equiripple) to the tangent, we can see that the 
resultant phase approximation is not exactly equiripple. 1 " We have not 



* Each design only required a few (e.g., 5) iterations. 

* We can see from Fig. 2 that the effect that the tangent transformation has on 
the error curve also depends on the values of the desired function. 
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implemented an algorithm to find the minimax solution to the phase, 
since, for our needs, the improvement in the phase approximation from 
the method outlined here did not seem to justify the use of a modified 
algorithm. 

III. A GRADIENT SEARCH TECHNIQUE FOR LEAST-SQUARES DESIGN 

The next design algorithm we describe involves the computation of 
the gradient vector relative to the set of coefficients in a product of 
quadratic factors. The transfer function of an all-pass digital filter, 
expressed as a product of second-order sections, is : 

The least-squares form 

E = E IPUk) - AjagH(ff-^f-)Jw(f k ) (12) 

will be used as a measure of the approximation error from the desired 
function D(f) on the set of frequency samples {/*}r\ Here, w(f) 
denotes a nonnegative weighting function. A. G. Deczky has also 
considered gradient techniques applied to the least-square design of 
all-pass digital filters. 7 In that paper, the emphasis was on envelope 
delay design. However, as shown in Section I, there are applications 
where envelope delay approximations are not adequate. Specifically, 
there are cases where phase distortion (e.g., constant phase offset) 
must be eliminated with an all-pass structure. 

Our design algorithm stems from familiarity with Ref. 8, which 
considers magnitude-only designs. With the least-squares criterion, 
the cascade second-order section form can be used. The advantage 
is that coefficient accuracy problems are minimized. As an alternative 
to the linear programming approach considered in Section II, this 
least-squares approach also enables one to more easily control the 
linear-phase offset permitted in the design. However, a disadvantage 
of the least-squares approach is that stability of the designed filter 
cannot be handled directly. Stability is obtained by confining the 
gradient movement to within the unit circle. This constraint may 
increase the likelihood of reaching an unsatisfactory local optimum. As 
we see later, there are initial guess procedures that provide excellent 
approximations to the desired phase function which, through the design 
algorithm, increase the likelihood of reaching a satisfactory local 
optimum. 
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3.1 Gradient calculations 

We find the entries of the gradient vector are 

f = - 2 t CD(A) - An gH (*—)>(/.) iAngfl^ 

da, * = i da, 

£--«£: COW - Ang ^.JXA) aAng f fl (C " Wt) • (14) 

op,- * = i op,- 

Here we define <*>(/) = Ang H(e~^) = tan" 1 / (/)/#(/), where 
/(/) = Im&gH(e-*") and R(f) = Real H(fi-**t). We seek 

^ - «(/)/«(/) - i(f)RUf) 

^f = RU)I,U) - l(f)RHf), 

where prime (') denotes the partial derivative relative to the subscript. 
After some algebra, we find 

|£ - 2(1 - 0i)Fi(f) sin 2*-/ i^i^M (15) 

OOti 

^ = 2ZW) (sin 4*-/ + a,- sin 2rf) 1 ^ i ^ M, (16) 

op, 

where F,(/) = |1 + *&-*** + $4T*"\-*. Finally, (13) and (14) can 
be simplified for 1 2j> i s£ M to 

|^ = - 4(1 - ft) E [D(/ t ) - *(/*)II*V(/*M/*) sin 2*/» (17) 

da,- it 

|| = - 4 E [!>(/*) - 0(A) iFiifowiU) (sin 4*/* + a, sin 2ift). (18) 
dpi k 

The minimization of E in (12) then proceeds with an iterative algorithm 
that is based on the formula 

C (n) = C (n-1> _ e B _ 1 il ll (V^) B _ lj (19) 

where c (n) is the coefficient vector (oi, (Si, a 2 , /3 2 , • • •, ct M , 0m) at the 
nth iteration, c„ is the nth step size in the coefficient adjustment, A n 
is a positive definite matrix at the nth step ( = / in the case of the 
steepest descent algorithm) and (VE)„ is the gradient vector whose 
entries are given by (17), (18) at the nth iteration (we use the Fletcher- 
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Powell algorithm). An initial guess procedure is required to start an 
iterative algorithm such as that of (19). 

3.2 Initial guess procedures for all-pass networks 

Convergence to a local minimum at which the approximation to a 
desired phase function is satisfactory can be made easier if a good 
initial guess is provided to c (0) of (19). A desired feature of an initial 
guess procedure is that it be simple in nature. After all, excessive 
computation and effort should not be expected in simply starting a 
complex algorithm. In this section, we consider two procedures in 
which only a linear set of equations need be solved to obtain initial 
values for \bk}k=o of (1). The value of having several initial guess 
procedures is that the designer may want to exercise his algorithm 
from multiple starting points to choose the best from a set of local 
optima. The following initial guess procedures operate on the direct 
form of i/(z -1 )(l) which can be factored to the cascade form (11). 

3.2.7 Tangent approximation by Gauss' method 

From (4) in Section II we know that a desired phase function can 
be approximated by considering a monotonic function of the phase, 
namely the tangent. Hence, 

N 

E bk sin 2irkj 

tan *(/) = " "-W (20) 

E b k cos 2irkf 

is the approximating function of the tangent of half the desired phase. 
If we require the estimates of the desired phase tangent [tan #<*(/)] to 
be "good" at a number of frequencies, we then have the following 
equations : 

AT N 

tan d (/o) E & fc cos 2tt/o — E b k sin 2irkf = r 

k=Q fc-1 

N N 

tan <t> d (fi) E &k cos 2irkfi - E b k sin 2irkfi = n 

k=Q k=l 



N N 

tan </>d(/ A /) E b k cos 2irkf M - E &* sin 2irkf M = r M . (21) 

/:=0 k-1 

If \r n }o were all zero, then the approximation would be exact. The 
objective then is to minimize E^=o rl, where M > N. This problem 
is a least-squares minimization problem for which the solution is 
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derived from solving a set of normal equations : 



or 



r (ai, ai)(ai, a 2 ) ■ 


• (ai, a*)] 


f&i] 




\(*i,dV 


(a 2 , ai)(a 2 , a 2 ) • 


■ ■ (a 2 , &n) 


b 2 


= 


(a 2 , d) 


I (a*, ai) 


(a. N , a. N ) 
Ah = e, 


V. 




.(a^d)^ 



(22) 



where 

a„ = (tan faifo) cos 2wnf — sin 2irnf , 

•tan (f>d(fi) cos 2irnfi — sin 2irnf h • • •, 

■tan 0d(/A/) cos 2-kuJm — sin 2Tnf M ) n = 1,2, 



N 



and 



Let 



d = [- tan d (/o), — tan 4> d (fi), • • •, - tan d (/^)] 
b = (61, 6 2 , • • •, hf) and e = [(ai, d), • • •, (a*, d)]. 



Pmax 



max 



\rn\ 



and 






where r* = (/"o, rt, • • •, rj/), the residual values after the least-squares 
approximation. If p max — p is large (it is always positive), then a 

Chebyshev approximation may be desirable. 18 

3.2.2 Tangent approximation in Chebyshev sense 

It is well known that the minimax solution to (21) requires solving 
an appropriate subsystem of N + 1 equations. Further, the minimax 
solution of N + 1 inconsistent equations can be effected by examining 
the least-squares solution to the same set of equations and proceeding 
to solve a set of N linear equations. 19 

For our purposes here, an effective method of obtaining an initial 
guess for the iterative procedure implied by (19) is that of choosing 
M = N and obtaining the minimax solution to (21). This can be done 
by solving (22) for b = (61, 6 2 , • • • , 6#) and then evaluating (21) for 
the residuals rl, r\, • • • , r* N . The minimax solution to (21) is then given 
by the linear set of equations 



Bb = 



(23) 



where B = (bj k ), N + 1 by N matrix with b jk = tan 4>d(fj) cos 2vkfj 
— sin 2-irkjj, a = e[sign (r ), sign (n), • • ■, sign (?•#)], and 



N N 

e=zZr?/Z\rt\. 

k=0 k—0 
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It may be noted that only N of N + 1 equations are used in the solu- 
tion of (23). 

3.2.3 Discussion 

It should be noted that no constraint has been made on the initial 
guess procedures of A or B to ensure that the resulting digital filter 
is stable. In fact, if £&* &* cos k2vf should ever change sign in |/| ^ \ 
or at least in the subinterval of approximation [/ ,/m], tnen tne 
transition from (20) to (21) is not really valid since a division by zero 
is implied. Should LtU b k cos k2rf be strictly positive over |/| ^ §, 
then stability results. 5 (The interesting point is that stability can 
result even if the cosine series does change sign in |/| ^ ?). However, 
the point to remember is that the resulting initial guess may be 
unstable. In our experience, we have not encountered any serious 
problems using these initial guess procedures. 

We must further remark that the inherent N sample delay present 
in these approximations [see (2)] could present a problem when 
designing filters with M ^ N sample delays. However, we feel, 
intuitively, that since some delay is unavoidable, a delay of the order 
of the filter will not, for most applications, be overly restrictive. 

The last point to consider is that the initial guess procedure of Sec- 
tions 3.2.1 and 3.2.2 obtains a direct form estimate of the digital filter 
coefficients. What is really required for c (0) of (19) are quadratic factors. 
We remark that we make the transition from the direct form estimate 
of (20) to quadratic factors by using a Bairstow quadratic factoriza- 
tion routine. 

3.3 Some considerations for least-squares design 

Often the engineering systems requirement of a digital filter can 
tolerate a linear-phase offset. While the systems engineer cannot 
always adapt to an arbitrary delay, there will usually be a range of 
delays permissible to him. How then can a designer incorporate these 
relaxations into the design mechanism? One technique for doing this 
is to add an acceptable delay to the desired function to create a new 
desired function and proceed from there. By designing filters for each 
of the permissible delays, one can choose from among the delays and 
their associated errors to decide which filter to implement. 

In Fig. 5 we can see the error function of a sixth-order filter* (B) 



* We have not tested the limit of the order of filters that can be designed by this 
method, but we have obtained a twentieth-order approximation (20-sample delay) 
to the desired function in this example. Quality initial guess procedures help us do 
this without excessive computation times. 
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Fig. 5 — Error curves for initial and final sixth-order Hilbert transformer designs. 



designed for a delay of six samples. The desired function is a 90-degree 
phase shifting filter with the approximation having weight 1 on 
[0.08, 0.41] and otherwise. Note the quality of the initial approxi- 
mation (A) using the first initial guess technique outlined in Section 
3.2. Of course, the disadvantage of presetting the delay is obvious; 
the choice of the optimal delay from those that are acceptable is not 
automatic but requires a separate design for each delay. However, 
eq. (12) can be expanded to include delay as parameter A 

E = E LD(fk) ~ AngH(e-*"») - A2rf k Jw(f k ). 

k 

An optimal A can be found analytically at each step in the gradient 
search and at convergence A will represent the amount of delay which, 
in conjunction with the filter, produces the best design.* Of course, we 
cannot expect that this delay will represent an integral number of 
samples or even a delay that the designer can tolerate. Figure 6 shows 
a desired function (A) (this curve is only shown where the weight of 
the approximation is nonzero) and its fourth-order approximation (B), 



* It is possible to include a constant phase angle as parameter B similar to the A 
used here. In such a case, our procedure becomes an envelope delay design technique. 
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Fig. 6 — Fourth-order approximation to desired shape (A). 



which is by solving for optimal A. Allowing for arbitrary delay, the 
algorithm obtained this optimal design with a 4.1-sample delay. 

We can offer an heuristic solution to guarantee integer delays in an 
automatic fashion ; namely, at each step (that is, at each calculation of 
A), the nearest acceptable delay* is used to replace A in the algorithm. 
This, of course, places a serious strain on optimality, although it does 
permit an automatic design procedure. 

As a footnote to this algorithm, we remark that there is a tendency, 
when working with procedures for designing filters in the cascade 
form, to use a previous optimal design of order n as the initial starting 
point in the design of filters of order n + 2. In the case of magnitude- 
only design, this is easily implemented since the appended second-order 
section can be initialized with magnitude 1. However, in the all-pass 
presentation there does not exist any second-order section that can be 
added which does not distort the overall phase when using a previous 
optimal design of order n to provide the initial guess for a design of 
order n + 2. And so the user of this algorithm must consider the effect 
of the appended second-order section if he does not want to obviate 
the value of a previous design toward providing an initial guess. 



* Nearest in 
"integer." 



the sense of greatest reduction of (12); "acceptable" here means 
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